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Abstract: In this note we develop an extension of the Marcenko-Pastur the- 
orem to time series model with temporal correlations. The limiting spectral 
distribution (LSD) of the sample covariance matrix is characterised by an 
explicit equation for its Stieltjes transform depending on the spectral den- 
sity of the time series. A numerical algorithm is then given to compute the 
density functions of these LSD's. 
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1. Introduction 

Let {X.j},j = 1, . . . , n be a sequence of p-dimensional real- valued random vectors 
and consider the associated empirical covariance matrix 



The study of the empirical spectral distribution (ESD) F n of S n , i.e. the distri- 
bution generated by its (real-valued) eigenvalues, goes back to Wishart in 1920's. 
A milestone work by Marcenko and Pastur (1967) states that if both sample size 
n and data dimension p proportionally grow to infinity such that \\mp/n = c for 
some positive c > and all the coordinates of all the vectors Xj's are i.i.d. with 
mean zero and variance 1, then F n converges to a nonrandom distribution. This 
limiting spectral distribution (LSD), named after them as the Marcenko-Pastur 
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distribution of index c has a density function 

fix) = -?—y/(b-x)(x-a), a<x<b, (2) 

ZTTCX 

with a = (1 — -y/c) 2 and b = (l + y^) 2 defining the support interval and has a point 
mass 1 — 1/c at the origin if c > 1 . Further refinements are made successively by 
many researchers including Jonsson (1982), Wachter (1978) and Yin (1986). 

An important work by Silverstein (1995) aimed at relaxing the independence 
structure between the coordinates of the Xj's and considered random vectors of 
form Yj = Tp^Xj where (T p ) is a sequence of non-negative definite matrices. 
Assuming that (T p ) is bounded in spectral norm and the sequence of ESD of 
(T p ) has a weak limit H, he established a (strong) LSD for the sample covari- 
ance matrix n" 1 ^YjYj and provides a characteristic equation for its Stieltjes 
transform. Despite a big step made by this generalisation, it still does not cover 
all possible correlation patterns of coordinates. Pursuing these efforts, a recent 
work by Bai and Zhou (2008) pushes a step further Silverstein's result by allow- 
ing a very general pattern for correlations between the coordinates of the Xj's 
satisfying a mild moment conditions. 

In this work, we extend such Marcenko-Pastur type theorems along another 
direction by considering time series observations instead of an i.i.d. sample. Let 
us first consider an univariate real-valued linear process 

oo 

zt = ^(j)k£t-k, teZ, (3) 

k=0 

where (e^) is a real- valued and weakly stationary white noise with mean zero and 
variance 1. The p-dimensional process (X() considered in this paper will be made 
by p independent copies of the linear process (zt), i.e. for Xj = [X\t, ■ ■ ■ , X p t) T , 

oo 

x it = ^2 0k£i,t-k, t e z, (4) 

fc=0 

where the p coordinate processes {(ei,t, • • • ,£p,t)} are independent copies of the 
univariate error process {et} in (3). Let Xi, . . . , X n be the observations of the 
time series at time epochs t = 1, . . . ,n. Again we are interested in the ESD of 
the sample covariance matrix S n in (1). 

The author should mention that a similar problem has been considered in Jin 
et al. (2009). However we propose much more general results in this note since 
firstly their results are limited to ARMA-type processes instead of a general 
linear process considered here and secondly, they do not find a general equation 
as the one proposed in Theorem 1 below except for two simplest particular cases 
of AR(1) and MA(1). 
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2. A Marcenko-Pastur type theorem for linear processes 

Recall that the Stieltjes transform of a probability measure \i on the real line 
is a map from the set C + of complex numbers with positive imaginary part onto 
itself and defined by 

Sn(z) = [ fj,(dx), z G C + . 

J x- z 

We always employ an usual convention that for any complex number z, \fz 
denotes its square root with nonnegative imaginary part. 

Theorem 1. Assume that the following conditions hold: 

1. The dimensions p — > oo, n — > oo andp/n — > c G (0,oo); 

2. The error process has a fourth moment: Kef < oo; 

oo 

3. The linear filter {<t>k) is absolutely summable, i.e. |<ftfc| < oo. 

fc=o 

Then almost surely the ESD of S n tends to a non-random probability distribution 
F. Moreover, the Stieltjes transform s = s{z) of F (as a mapping from C + into 
C + ) satisfies the equation 

1 1 f 2w 1 

z = 1 / ? — r d\ , (5) 

where /(A) is the spectral density of the linear process (zt): 



> 



k=0 



, AG[0,2vr). (6) 



The proof of the theorem is postponed to Section 4. Let us mention that 
although the case c = is beyond the scope of Theorem 1, Equation (5) leads in 
this case to the solution s(z) = 1/(70 — z), that is the LSD would be the Dirac 
mass at 70 = Var(X#). This conjectures an extension of Theorem 1 to the so- 
called "very large p and small n" asymptotics where one assumes p — > 00, n — > 00 
and p/n — > 0. Indeed, in this scenario taking into account that the population 
covariance matrix of Xj equals 70 ip, one can expect that the sample eigenvalues 
of S n stay close to the population ones (all equal to 70). Note that such results 
exist for i.i.d. sequence (Xf) with i.i.d. components (see Bai and Yin, 1988). 
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Figure 1 . Determination of the support [x\ , asa] of the LSD. Left panel c > 1 and right panel 
c < 1. 



2.1. Support of the LSD F 

Starting from Eq. 5 and following the techniques devised in Silverstein and Choi 
(1995), we can describe precisely the support of the LSD F in previous theorem. 

Let a and b be respectively the minimum and the maximum of the function 
2nf over [0, 2tt] . As / is infinitely differentiable and positive everywhere, both 
a and b are attained and the range of 2irf is exactly the interval [a, b]. We 
will always exclude the situation a = since it corresponds to a special class 
of linear processes, namely non-invertible ARMA processes, see Grenander and 
Szego (1958, chap. 9), which has no practical interest for applications. Therefore 
the map s i— > z = g(s) in Eq.(5) has a trace for real-valued s providing s ^ 
[— — , — r-]. Figure 1 depicts this map for both c < 1 and c > 1 cases. 

The following proposition is a straightforward application of results from 
Silverstein and Choi (1995) and we then omit its proof. 

Proposition 1. With the map g : s \— > z = g(s) in Eq.(5) restricted to real 
s & I - ac> — ?H (figure 1) the following holds: 

1. The LSD F has a compact support \x\,X2[ C [0,oo) on which it has a 
continuous density function. In case of c > 1, F has an additional point 
mass 1 — 1/ c at the origin. 

2. When c > 1, the map g has an unique maximum s\ on (— oo, — -M and 
an unique minimum S2 on (— ^,0) and we have The edges of the support 
interval are given by these local extrema: x\ = g(s\) and X2 = g{s2)- 

3. When c < 1, the map g has an unique maximum s\ on (0, oo) and an 
unique minimum S2 on (— ip,0). The edges of the support interval (for the 
absolutely continuous component) are again given by these local extrema: 
x\ = g(si) and x 2 = g(s 2 ). 
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2.2. Application to an ARM A process 

For simplicity, we consider the simplest causal ARMA(1,1) process for the coor- 
dinates: 

zt = (pzt-i + e t + 9e t -i, t £ Z, 

where \(j>\ < 1 and 9 is real. The aim is to find a simplified form of general 
equation (5). We have 

1 l-6e iXl2 



2tt/(A) 

and 



1 + 0e % 



If 2 * 1 1 / 1 d£ 

dX - 



2vr J cs + {2vr/(A)}- 1 2vri f m=1 cg + 



2 



By a lengthy but elementary calculation of residues detailed in Section 4, we find 

l vy + g)(l + ^) e (a) \ 



cse-4>\ 6{cs9 - 4>) J ' 

with 

cs (l + 9 2 ) + 1 + (j> 2 
a= cs e-(j) ' e(a) = sgn(9a) . (8) 

Therefore for an ARM A (1,1) process, the general equation (5) reduces to 

= _\ _J_ (0 + 0)(l + 00) 6(a) 

a cs9 - ( CS - (p) 2 ■ W 

Let us mention that it is important to have an explicit formula for the integral 
in (5) to implement numerical algorithms like the one proposed in Section 3 in 
order to compute the density function of the LSD F. 

Case of an AR(1). For this particular case, we have 9 = and a = —(cs + 
1 + 4> 2 )/<i). As 9s > 0, e(a) = sgn(-0). It follows that 

(<i> + 9)(l + (t>9) e(a) sgn(-0) 1 



(cs9-(p) 2 Va 2 -4 -<t> / (cs+i+W ~i 

1 1 1 

W\ / (cs+i+^r 4 ~ ^(cs + 1 + 2 ) 2 - 40 2 ' 



Therefore the Stieltjes transform s of the LSD is solution to a simpler equation 

z = — + . . (10 

s ^(cs + 1 + 02)2 _ 402 
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It is worth noticing that if we further assume <j> = 0, this equation reduces to 
z = —1/s + l/(cs + 1) which characterises the standard Marcenko-Pastur law 
with i.i.d. coordinates. Furthermore, for the determination of the support [xi,X2] 
of the LSD, we notice that 

so that its extrema are a = 1/(1 + \4>\) 2 and b = 1/(1 — \4>\) 2 (see Figure 1). 



Case of an MA(1). Here we have 4> = and 



e \cs 

Hence e(a) = — sgn(#) and it is readily checked out that the Stieltjes transform 
of the LSD is solution to the equation 

111 1 

z = — + — + ^T2 i ( n ) 

S CS CS J(l + 1 + 02 f_ 4d 2 

Again if we further assume 9 = 0, this equation reduces to the one for the stan- 
dard Marcenko-Pastur law. Furthermore, for the determination of the support 
[xi, Xi\ of the LSD, we notice that 

2vr/(A) = |l + 6>e iA | 2 , 

so that its extrema are a = (1 — \9\) 2 and b = (1 + \9\) 2 (see Figure 1). 



3. A numerical method for computing the LSD density function 

In this section we provide a numerical algorithm for the computation of the 
density function h of the LSD defined in Eq.(5) through its Stieltjes transform 
s. We have 

1 

S ~ -z + A(s{z)) 

with 

1 f 27T 1 

= 7T / ^ J9 fmw dX . 
2tt J cs + {2tt/(A)} 1 

The algorithm we propose is of fixed-point type. 
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Algorithm For a given real x, let e be small enough positive value and set 
z = x + is. 

Choose an initial value sq{z) = u + ie and iterate for k > the above 
mapping 

Sk+i(z) = {-z + A(s k {z))Y 1 , 

until convergence and let sk{z) be the final value. 
Define the estimate of the density function h(x) to be 

h(x) = —Qsk(z) . ■ 

It is well-known that this iterated map has good contraction properties guar- 
anteeing the convergence of the algorithm. There are however two issues which 
need a careful consideration. First the integral operator A is usually approxi- 
mated by a numeric routine and because of a high number of calls to A, the 
resulting algorithm is slow. In this aspect analytic formula for A when available 
are well acknowledged as Eq.(9) in the case of an ARMA(1,1). 

A second issue is that overall we first need to determine the support interval 
[xi, x<i\ of the density function h. This is handled with the help of description of 
x/s given in Proposition 1. 

For four ARM A (1,1) models listed in Table 1, we have used this algorithm 
with the map A defined in (9) to get the density plots displayed in Figure 2. 

Table 1 

ARM A( 1,1) models for density plots. Reference support interval for Marcenko-Pastur law of 



index c 


= 0.2 is [0.306,2.094]. 


Parameters (<f>, 8, 


c) Estimated support [asijXa] 


(0.4, 0, 0.2) 


[0.310, 2.875] 


(0.4, 0.2, 0.2) 


[0.319, 3.737] 


(0.4, 0.6, 0.2) 


[0.382, 6.186] 


(0.8, 0.2, 0.2) 


[0.485, 13.66] 



Compared to the reference standard Marcenko-Pastur law with the same 
dimension to sample ratio c, the above density functions from ARMA models 
share a similar shape with however a support interval getting larger and larger 
with increasing ARMA coefficients <f> and 6. 
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4. Proofs 

Proof of Theorem 1 

Recall that the p coordinates of the vectors Xi, . . . , X n are i.i.d. while the tem- 
poral covariances Cov(Xj s , Xn) are by definition those of (zt)' for all 1 < i < p, 

Cov(X is , X it ) = Cov(z s ,z t ) = 7t_ s , 1 < s,t < n . 

Let / = 27r/. It follows that the covariance matrix T n of each coordinate process 
(Xn., . . . , Xi n ) equals to the n-th order Toeplitz matrix associated to /: 

T n (s,t) = 7t_ s , l<s,t<n, 

and 

oo 

/(A) = Y, ^e tkX , Ag[0,2tt). 

k=~ oo 

We are going to apply Theorem 1.1 of Bai and Zhou (2008) for a strong 

n 

limit of the ESD of the sample covariance matrix S n = jjr Yl XtXj . Under the 

t=i 

assumptions made, all the conditions of this theorem are satisfied except that we 
need to ensure a weak limit for spectral distributions of (T n ). 

The function / belongs to the Wiener class, i.e. the sequence of its Fourier 
coefficients is absolutely summable. Moreover note that / is infinitely differen- 
tiable, its minimum a and maximum b are attained. According to the fundamental 
eigenvalue distribution theorem of Szego for Toeplitz forms, see Grenander and 
Szego (1958, sect. 5.2), for any function (p continuous on [a, b] and denoting the 

(n) (n) 

eigenvalues of T n by a{ , . . . , &n , it holds that 

hm -X>(af) = 1- f\{J(\))d\ . 

Consequently, the ESD of T n (i.e. distribution generated by the ai^'s) weakly 
converges to a nonrandom distribution H with support [a, b] and defined by 

H(x) = — [ x d\ , (12) 

and we have for ip as above, 

^ ^x)dH(x) = — ip(f(X))d\ . (13) 
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Furthermore, by application of Theorem 1.1 of Bai and Zhou (2008), it holds 
that the ESD of ^S n converges almost surely to a nonrandom probability distri- 
bution whose Stieltjes transform m solves the equation 

z = -- + - I -^—dH{x) 
m c J 1 + mx 

= + =— / , 

m 2-kcJq m + l/f 

where we have used (13) in the last equation. The equation(5) follows by observ- 
ing the relation s(z) = ^m(z/c). 

Proof of Equation (7) 

The aim is to evaluate the integral 

If 1 



1 2-/ 



l«l=i 



cs + 



1+0$ 



Let m = cs in this computation of residues. We have 

i r l + ^ + ^ + r 1 ) <% 



2ni J m=l m {1 + 9 2 + 0(f + £-*)} + {1 + <P 2 ~ <K£ + f" 1 )} £ 

i / fi (0 + e)(i + ^) i 



m9-<j)27ri J\ &=1 \£ 9{m9 - <j>) £ 2 + 1 + a£ 
> + 0)(l + 00) 1 / d£ 1 



1 



with 



m0-<j)\ 9{m9 - (j>) 2iri J\^ =1 P(£) j ' 



n/^ .2 . 1 m(l + # 2 ) + 1 + </> 2 



m# — 4> 

Let £ x , £ 2 be the roots of P(£) = £ 2 + 1 + af . Then 

As ^1^2 = 1, only one of the two poles is inside the unit circle. It is readily checked 
that if 9a > 0, then < |^| and 

I £ d£ 1 1 



2tu J n=1 P(£) 6-6 Vc^4' 
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Otherwise we have |£i| > |^| and the integral has an opposite sign. Summarising 
both cases we get 

1 r d£ e(a) 
with e(a) = sgn(9(a)). Equation (7) is proved. 
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References 

Bai, Z., Zhou, W., 2008. Large sample covariance matrices without independence 

structures in columns. Statist. Sinica 18 (2), 425-442. 
Bai, Z. D., Yin, Y. Q., 1988. A convergence to the semicircle law. Ann. Probab. 

16 (2), 863-875. 

Grenander, U., Szego, G., 1958. Toeplitz forms and their applications. California 
Monographs in Mathematical Sciences. University of California Press, Berkeley. 

Jin, B., Wang, C, Miao, B., Lo Huang, M., 2009. Limiting spectral distribu- 
tion of large-dimensional sample covariance matrices generated by varma. J. 
Multivariate Anal. 100, 2112-2125. 

Jonsson, D., 1982. Some limit theorems for the eigenvalues of a sample covariance 
matrix. J. Multivariate Anal. 12 (1), 1-38. 

Marcenko, V., Pastur, L., 1967. Distribution of eigenvalues for some sets of ran- 
dom matrices. Math. USSR-Sb 1, 457-483. 

Silverstein, J. W., 1995. Strong convergence of the empirical distribution of eigen- 
values of large-dimensional random matrices. J. Multivariate Anal. 55 (2), 331- 
339. 

Silverstein, J. W., Choi, S.-L, 1995. Analysis of the limiting spectral distribution 
of large-dimensional random matrices. J. Multivariate Anal. 54 (2), 295-309. 

Wachter, K. W., 1978. The strong limits of random matrix spectra for sample 
matrices of independent elements. Ann. Probability 6 (1), 1-18. 

Yin, Y. Q., 1986. Limiting spectral distribution for a class of random matrices. 
J. Multivariate Anal. 20 (1), 50-68. 



imsart -generic ver. 2009/08/13 file: lp-ss.tex date: September 9, 2011 



J. Yao/Marcenko-Pastur theorem for time series 



11 




imsart -generic ver. 2009/08/13 file: lp-ss.tex date: September 9, 2011 



